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Abstract We derive an analytical approximation of nonlinear force-free mag- 
netic field solutions (NLFFF) that can efficiently be used for fast forward- 
fitting to solar magnetic data, constrained either by observed linc-of-siglit mag- 
netograms and stereoscopically triangulated coronal loops, or by 3D vector- 
magnetograph data. The derived NLFFF solutions provide the magnetic field 
components Bx{:x.), By{x.), i3^(x), the force- free parameter a(x), the electric 
current density j(x), and arc accurate to second-order (of the nonlinear force- 
free a-parametcr) . The explicit expressions of a force- free field can easily be 
applied to modeling or forward-fitting of many coronal phenomena. 
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1. Introduction 

The coronal magnetic field can be constrained in a number of ways, such as 
by extrapolation of photospheric magnetograms or vector-magnetograph data, 
by radio observations of gyroresonance layers above sunspots, of by coronal 
seismology of oscillating loops. Before the advent of the STEREO mission, 
attempts were made to model observed coronal loops with stretched potential 
field solutions (Gary and Alexander, 1999), to fit a linear force- free model with 
solar-rotation stereoscopy (Wiegelmann and Ncukirch, 2002; Feng et ai, 2007a); 
by tomographic reconstruction with magnetohydrostatic constraints (Wiegel- 
mann and Inhester, 2003; Ruan et ai, 2008), by magnetic modeling applied to 
spectropolarimetric loop detections (Wiegelmann et ai, 2005), or by magnetic 
field supported stereoscopic loop triangulation (Wiegelmann and Inhester, 2006; 
Conlon and Gallagher, 2010). Recently, stereoscopic triangulation of coronal 
loops with the STEREO mission became available, which constrains the 3D 
geometry of coronal magnetic field lines (Aschwanden et al., 2008; Aschwanden, 
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2009). The plethora of coronal high- resolution data allows us now to compare 
different magnetic models and to test whether they are self-consistent. A critical 
assessment of nonlinear force-free field (NLFFF) codes revealed the disturbing 
fact that different NLFFF codes yield incompatible results among themselves, 
and exhibit significant misalignments with stereoscopically triangulated loops 
(DeRosa et ai, 2009; Sandman et ai, 2009; Aschwanden and Sandman, 2010; 
Sandman and Aschwanden, 2010; Aschwanden et ai, 2012a, b). The discrep- 
ancy was attributed to uncertainties in the boundary conditions as well as 
to the non-forcefreeness of the photosphere and lower chromosphere. Earlier 
tests with the virial theorem already indicated that the magnetic fields in the 
lower chromosphere at altitudes of ft. < 400 km are not force-free (Metcalf et 
al., 1996). Constraints by coronal tracers thus have become an important cri- 
terion to bootstrap a self-consistent magnetic field solution. The misalignment 
between theoretical extrapolation models and stereoscopically triangulated loops 
could be minimized by using potential field models with forward-fitted unipolar 
magnetic charges (Aschwanden and Sandman, 2010) or dipoles (Sandman and 
Aschwanden, 2011). 

In this Paper we go a step further by deriving a simple analytical approxima- 
tion of nonlinear force-free field solutions that is suitable for fast forward-fitting 
to stereoscopically triangulated loops or to some other coronal observations. 
While accurate solutions of force- free magnetic fields have been known for special 
mathematical functions (Low and Lou, 1990) that have been used to reconstruct 
the local twist of coronal loops (Malanushenko et al., 2009, 2011), they are not 
suitable for forward-fitting to entire active regions. In contrast, our theoretical 
framework entails the representation of a potential or non-potential field by a 
superposition of a finite number of elementary field components that are associ- 
ated with buried unipolar magnetic charges at arbitrary locations, each one being 
divergence-free and force-free to a good approximation, as we test numerically. 
While this Paper contains the analytical framework of the magnetic field model, 
the numerical forward-fitting code with applications to observations will be pre- 
sented in a Paper II (Aschwanden and Malanushenko, 2012), and applications 
to stereoscopically observed active regions in Aschwanden et al, (2012a, b). 



2. Theory 

2.1. Potential Field Parameterization 

The simplest representation of a magnetic potential field that fulfills Maxwell's 
divergence-free condition (V • B = 0) is a unipolar magnetic charge j that is 
buried below the solar surface, which predicts a magnetic field B j (x) that points 
away from the buried unipolar charge and whose field strength falls off with the 
square of the distance rj , 




(1) 
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where Bj is the magnetic field strength at the solar surface above a buried 
magnetic charge, {xj,yj, Zj) is the subphotospheric position of the buried charge, 
dj is the depth of the magnetic charge, 

d^ = ^-^^' + y' + -h (2) 

and Vj ~ [x ~ Xj,y — yj, z — Zj] is the vector between an arbitrary location x = 
{x, y, z) in the solar corona (were we desire to calculate the magnetic field) and 
the location {xj,yj,Zj) of the buried charge. We choose a Cartesian coordinate 
system (x, y, z) with the origin in the Sun center and are using units of solar 
radii, with the direction of z chosen along the line-of-sight from Earth to Sun 
center. For a location near disk center (x <C 1,?/ <C 1), the magnetic charge 
depth is dj « (1 — Zj). Thus, the distance rj from the magnetic charge is 

rj ^ \J{x- Xjf + {y- yjY + {z - z^Y . (3) 

The absolute value of the magnetic field Bj {rj ) is simply a function of the radial 
distance rj (with Bj and dj being constants for a given magnetic charge), 

S('-.) = S.(^)' ■ (4) 

In order to obtain the Cartesian coordinates {B^, By, B^) of the magnetic 
field vector Bj(x), we can rewrite Equation (1) as, 

Bx{x,y,z) = Bj {dj/rj f {x - Xj)/rj 

Byix, y, z) = Bj [dj/rjf {y - yj)/rj . (5) 
Bz{x, y, z) = Bj {dj/rjf {z ~ Zj)/rj 

We progress now from a single magnetic charge to an arbitrary number iVm of 
magnetic charges and represent the general magnetic field with a superposition 
of N-ai buried magnetic charges, so that the potential field can be represented by 
the superposition of N^^ fields Bj from each magnetic charge j = 1, ...,N^, 

b(x)=5:b,(x)=5:sj^ ^. (6) 

As an example we show the representation of a dipole with two magnetic 
unipolar charges (TVm = 2) of opposite polarity {B2 = —Bi) in Figure 1. 
Each of the unipolar charges has a radial magnetic field (dotted lines), but 
the superposition of the two vectors of both unipolar charges in every point of 
space, B(x) = Bi(x) + B2(x), reproduces the familiar dipole field. For the case 
shown in Figure 1 we used the parameterization of two subphotospheric unipolar 
magnetic charges at positions xi = —0.5 and X2 — +0.5, which produces dipole- 
like field lines (solid curves), while they converge to the classical solution of a 
dipole field in the limit of xi and X2 1— >• 0, as it can be shown analytically 
(Jackson, 1972; p.l84). 
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Figure 1. The magnetic field of a symmetric dipole (dashed fines) is shown, together with the 
field resulting from the superposition of two unipolar magnetic charges (solid lines). The two 
field models become identical once the two unipolar charges are moved towards the location 
of the dipole moment at position (x, y) = (0, 0). The radial field of each unipolar (positive and 
negative) charge is also shown for comparison (dotted lines). 



2.2. Force-Free Field Solution of a Uniformly Twisted Fluxtube 

A common geometrical concept is to characterize coronal loops with cylindrical 
fluxtubes. For thin fluxtubes, the curvature of coronal loops and the related 
forces can be neglected, so that a cylindrical geometry can be applied. Because 
the footpoints of coronal loops are anchored in the photosphere, where a random 
velocity field creates vortical motion on the coronal fluxtubes, they are gener- 
ally twisted. We consider now such twisted fluxtubes in a cylindrical geometry 
and derive a relation between the helical twist and the force-free parameter a. 
The analytical solution of a uniformly twisted flux tube is described in several 
textbooks {e.g., Gold and Hoylc, 1960; Priest, 1982; Sturrock, 1994; Boyd and 
Sanderson, 2003; Aschwanden, 2004), but we summarize the derivation here to 
provide physical insights for the generalized derivation of nonlinear force-free 
magnetic field solutions derived in Section 2.3 in a self-consistent notation. 

We consider a straight cylinder where a uniform twist is applied, so that an 
initially straight field line B = (0,0, Bs), aligned with a field line coordinate s, 
is rotated by a number A^twist of full turns over the cylinder length I, yielding an 
azimuthal field component at radius p, 

pd(f> 2npNt„ist , 

57 = ^ = — I — = 

with the constant b defined in terms of the number of full twisting turns iVtwist 
over a (loop) length I. 

The cylindrical geometry of a twisted flux tube is visualized in Figure 2. 
The longitudinal component of the untwisted magnetic field corresponds to a 
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Figure 2. The basic 3D geometry of a cylindrical flux tube with uniform twist is defined by 
the length I of the cylinder axis, the number of twisting turns along this length, Aftwisti or by 
the misalignment angle fj, at the flux tube radius p between the potential field line (aligned 
with the cylindrical axis) and the non-potential field line B'^^ (aligned with the twisted loop). 
The non-potential field line B'^^ can be decomposed into a longitudinal field component Bs 
and an azimuthal field component B^. 

potential field vector B^, while the twisted non-potential field line B^^ has 
a helical geometry with an angle /x at a radius p, which can be described by 
the longitudinal component Bg and the azimuthal component B^p. The fluxtube 
can be considered as a sequence of cylinders with radii p, each one twisted by 
the same twist angle dip/ds = 27riVtwist/'- For uniform twisting, the magnetic 
components B^p and Bs depend only on the radius p, but not on the length 
coordinate s or azimuth angle Lp. Thus, the functional dependence in cylindrical 
coordinates [p, Lp, s) is 



Consequently, the general expression of V x B in cylindrical coordinates. 



is simplified with Bp — Q and the sole dependencies of B^f, (p) and Bg (p) on the 
radius p (Equation (7)), yielding a force-free current density j of. 



Requiring that the Lorentz force is zero for a force-free solution, F = j x B = 0, 
we obtain a single non-zero component in the radial p-direction, since jp — 



B=[Bp,B^,Bs] = [0,B^{p),Bs{p)] . 



(8) 




(9) 




(10) 
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and Bp ~ for the two other components, 

F=j xB = [S,j^-B^j„0,0] , (11) 
yielding a single differential equation for Bs and B^, 

B,^+B^--f{pB^)^0. (12) 
dp pup 

By substituting B^p = bpBg from Equation (7) into Equation (12), this simplifies 
to, 

^[(l + 6V)Ss]=0. (13) 

A solution is found by making the expression inside the derivative to a constant 
(-Bo), which yields B^p and Bg, 



B — [Bp, Bp, Bs] 



Q Bo bp Bo 



1 + b^p^ ' 1 + 62 p2 



(14) 



[This equation corrects also a misprint in Equation (5.5.8) of Aschwandcn (2004; 
p. 216), where a superflous zero component has to be eliminated]. With the 
definition of the force-free a-parameter. 

An 

(V X B) = — j = a{p)B , (15) 
c 

we can now verify that the a-parameter for a uniformly twisted fluxtube depends 
only on the radius p, 

= (TW) ' ^'^^ 

with the constant b defined in terms of the number of full twisting turns A'twist 
over a (loop) length / (see Equation (7)), 

27rA'"twist .,„x 
6=—^. (17) 

The geometry of a twisted flux tube is visualized in Figure 3 (top middle), 
where the parallel field lines are aligned with the coordinate axis s in the vertical 
direction, the cross-sectional radius p is defined in the direction perpendicular 
to s, and the twist angle (p is indicated in the horizontal projection (Figure 
3, bottom middle). According to Equations (8) and (14), the variations of the 
longitudinal Bs{p) and of the azimuthal component Bp{p) with radius p are, 

Bsip) = ^ , (18) 



SOLA: ms.tex; 13 July 2012; 0:13; p. 6 



Nonlinear Force-Free Magnetic Field 





Figure 3. The field line geometry is shown for an untwisted cylindrical flux tube (left), a 
twisted cylindrical flux tube (middle), and for a twisted radial field (right), from the side view 
in the X2:- plane (top) and from the top view in the xj/- plane (bottom). The top panels show 
the longitudinal magnetic field component -Bs(p) and the bottom panels show the azimuthal 
magnetic field component By,(p,ip). 



These radial dependencies are shown in Figure 4 for different numbers of twist 
(-^twist = 0.5, 1.0, 1.5). In the limit of vanishing twist (A'^twist = i-^ 6 = 0), we 
have an untwisted flux tube (Figure 3 left) with a constant longitudinal field 
Bs{p) ~ Bq and a vanishing azimuthal component B^(p) = 0. The dependence 
of the azimuthal field component B^(p) and the longitudinal field component 
Bs{p) as a function of the radius p from the twist axis (Figure 4) shows that 
the longitudinal component falls off monotonically with radius p, while the az- 
imuthal component increases first for small distances p <^ I, but falls off at 
larger distances. Thus, the twisting causes a smaller cross-section of a fluxtube 
compared with the potential field situation, as widely known {e.g., Klimchuk et 
al, 2000). 



SOLA: ms.tex; 13 July 2012; 0:13; p. 7 



M.J. Aschwandcn 



at 




0.0 



0.0 



0.2 



0.4 



0.6 



0.8 



1.0 



radial distance p/l 



Figure 4. The dependence of the longitudinal (solid lines) and azimuthal magnetic field 
component (dashed lines) as a function of the distance r/l from the twist axis field is shown 
for three different amounts of twist (A^twist = 1-5, 1.0, 0.5 full turns per loop length I). 

2.3. Nonlinear Force-Free Field Parameterization 

We are now synthesizing the concept of point-like buried magnetic charges that 
we used to parameterize a potential field (Section 2.1) with the uniformly twisted 
flux tube concept that represents an exact solution of a nonlinear force-free field 
(Section 2.2). The geometric difference between the two concepts is the spherical 
symmetry of a point charge versus the parallel field configuration of an untwisted 
flux tube. However, we can synthesize the two geometries by considering the 
parallel field as a far-field approximation of a radial field. In an Euclidean parallel 
field, the equi-potential surface is a plane perpendicular to the parallel field 
vector, while a radial field has spherical equi-potential surface. We can make the 
transformation of a parallel field in cylindrical coordinates [s, p,Lp) into a radial 
field with spherical coordinates (r, 9, (p) by mapping (see Figure 3), 



This transformation from cylindrical to spherical coordinates preserves the or- 
thogonality of the longitudinal field component {Bs i-^- Br) to the equi-potential 
surface (s =const i— >■ r =const) and conserves the magnetic fiux $(?') along a 
bundle of field lines with area A{r) = p^ir), 



if the longitudinal component B{r) oc (Equation (1)) decreases quadratically 
with distance from the magnetic charge. Thus, applying the transformation 
into spherical coordinates (Equation (20)) and the magnetic flux conservation 
(Equation (21)) to the straight flux tube solution (Equations (18) and (19)), 



p H- r sin(^^) 



(20) 



$(r) B{r)A{r) = B{r)p^{r) = B{r)r'^ sin^ 9 = const , 



(21) 
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we can already guess the approximate nonlinear force-free solution in spherical 
coordinates, 



Brir,e) oc r- 



1 



(1 -I- 52 r 2 sin^ 



(22) 



br sin 6 



oc r 



(23) 



(1 + sin^ 61) ■ 

More rigorously, we can derive a nonlinear force-free field solution by writing 
the divergence-free condition (V • B) = and the force-free condition (V x B) = 
(47r/c)j = a(p)B (Equation (15)) of a magnetic field vector {Br, Bg, B^) in 
spherical coordinates (r, 6, ip) (with the origin at the location of the magnetic 
charge and the spherical symmetry axis aligned with the vertical direction to 
the local solar surface), 



(V.B) = l|-(r2i?. 



1 



d 



r sin 6 89 



[Bg sin 6*) 



1 dB^ 
r sin dip 



= 



(24) 



[V X B], 



1 



[V X B], = 



rsin^ 
1 



d 

^{B^ sin( 



1 dBr d 



sinO dip dr 



dB^ 

dip 



aBr 



aBg , 



(25) 
(26) 



fV X Bl,. = - 



a dBr 



= aB, 



(27) 



For a simple approximative nonlinear force- free solution we require axi-symmetry 
with no azimuthal dependence (djdip = 0) and neglect components that con- 
tribute only to second order (Bg oc [6rsin6']2 w 0), in analogy to the uniformly 
twisted flux tubes on cylindrical surfaces (Figure 2). This requirement simplifies 
Equations (24)-(27) to, 







1 



d 



r sin 9 d9 
1 d 



[B^p sin t 
(rB^) 



aBr 



r dr 
1 dB. 



r de 



yB, 



if ■ 



(28) 
(29) 
(30) 
(31) 



Eliminating a from Equations (29) and (31) and using the analog ansatz as for 
cylindrical fluxtubes (Equation (7)), 



B,r 



Brbr sin 9 



(32) 
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we obtain a similar differential equation as in Equation (13), 

^ [B,(l + 5Vsin2 0)] =0 . (33) 

A solution of this differential equation is obtained by setting the expression inside 
the bracket to the constant Bo{(P /r^), which fulfills the divergence- free condition 
(Equation (28)), and we obtain a solution for Br and (using Equation (32)), 
for a (using Equation (29)), 

Br{r, 9) = Bo f 4) ,^^,^\ . , (34) 
V? / (1 + o"^r"^ sm 6) 

, f <P\ br sin , , 

B^ir, e)=Boi-) ,^ , . . 35 

J (1 + b^r^ sm 0) 

Bg{r,e)^0, (36) 



26cos6> 
(1 + o"^r^ sui 0) 

This solution fulfills both the force- free condition (Equations (29)-(31)) and the 
divergence- free condition (Equation (28) ) to second-order accuracy (oc [br sin 6*] ^ ) .| 
We see that this solution is identical with the simplified derivation of Equations 
(22) and (23). At locations near the twist axis {9 i— > 0), the general solution 
(Equations (34)-(37)) converges to the cylindrical flux tube geometry solution 
(Equations (18) and (19)). Furthermore, in the limit of vanishing twist (6 n> 0) 
we retrieve the potential-field solution (Equation (4)), since the force-free pa- 
rameter becomes a n- 0, the azimuthal field component becomes B^ = 0, and 
the radial component reproduces the potential-field solution B,- M> Bo{(P/r^). 



2.4. Cartesian Coordinate Transformation 



In the derivation in the last section we derived the solution in terms of spherical 
coordinates (r, 9, ip) in a coordinate system where the rotational symmetry axis 
is aligned with the vertical to the solar surface intersecting a magnetic charge 
j. Since we are going to model a number of magnetic charges at arbitrary po- 
sitions on the solar disk, we have to transform an individual coordinate system 
{rj , 9j , tpj ) associated with magnetic charge j into a Cartesian coordinate sys- 
tem (x, y, z) that is given by the observers line-of-sight (in z-direction) and the 
observer's image coordinate system (x, y) in the planc-of-sky. The variables of 
the Cartesian coordinate transformation are shown in Figure 5. 

The radial magnetic field vector (which is pointing radially away from a 
magnetic charge j located in the solar interior at (xj,yj^Zj) is simply given by 
the difference of the Cartesian coordinates from an arbitrary location {x,y,z), 

= [cOSr,a:,COSr,y,COSr,2] , (38) 



By 

Br 



Xj y^yj z~Zj 
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Figure 5. The geometry of a twisted radial field of a magnetic charge j buried at a sub- 
photospheric position (xj ,yj , Zj) is shown. The central twist axis (dashed line) intersects an 
equi-potential surface at position {xo,yo, zo) and the longitudinal field vector Br at position 
[x, y, z) has a radial distance p from the twist axis and an azimuth angle ip. The azimuthal 
magnetic field component at location [x, y, z) is orthogonal to the radial vector p and the 
longitudinal field component Br, as well as to the direction of the twist axis R. 



where Br is the absolute value of the radial magnetic field component Br{rj, 0j) 
(Equation (34)), rj is the spatial length of the radial vector rj (Equation (3)), 
defining the directional cosines cos^,.; (for the 3D coordinates i — x,y, z) of the 
radial magnetic field vector B^. 

The azimuthal component Bj^ (with the absolute value B;p{rj,9j) defined in 
Equation (35)) of the twisted magnetic field is orthogonal to the direction of the 
twist axis R (aligned with the local vertical), 

R = [xj,yj,zj] , (39) 

and the radial magnetic field component B^ (Figure 5), and thus can be com- 
puted from the vector product of the two vectors B^ and R, 



Bj^ R X By^ 
B^ |R X B^ 



_ — ^ = [coSi^^a;,cos<p,j^,coSipJ , (40) 



which defines the directional cosines cos,^^i of the azimuthal component in the 
Cartesian coordinate system. The vector product allows us also to extract the 
inclination angle 9j between the radial magnetic field component B^ and the 
local vertical direction R, 

1 /|R X B^ 
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Finally, the total non-potential magnetic field vector B ~ [Bx-, By,Bz) is then 
the vector sum of the radial B,. and the azimuthal magnetic field component B;^, 

Ex = Br{rj,9j)coSr^x+B^{rj,6j)cos^,x 
By = Br{rj,6j)coSr,y+B^{rj,9j)cos^,y , (42) 
= Br{rj,9j) coSr,z +B^{r 6 j)cos^^z 

with the directional cosines (coSr,i, cos^^i, cose^i) defined by Equations (38) and 
(40). This is a convenient parameterization that allows us directly to calculate 
the magnetic field vector of the non-potential field Bj = {Bx, By, B^) asso- 
ciated with a magnetic charge j that is characterized with five parameters: 
{Bj , Xj , Uj , Zj ,aj), where we define the force- free a-parameter from the twist 
parameter bj = 27riVtwist/^ (Equation (7)) at the location of the twist axis 
{0, = 0), 

aj = a{Bj = 0) = 2 6j , (43) 

according to Equation (37). 

2.5. Superposition of Twisted Field Components 

The total non-potential magnetic field from all j — \, ...,N,^ magnetic charges 
can be approximately obtained from the vector sum of all components j (in an 
analog way as we applied in Equation (6) for the potential field), 

B(x) = gB,(x), (44) 

where the vector components Bj = (B^j, Byj, B^j) of the non-potential field 
of a magnetic charge j are defined in Equation (42), which can be parameterized 
with free parameters (Bj,Xj,yj,Zj,aj) for a non-potential field, or with 

4A^ni free parameters for a potential field (with aj =0). Of course, the sum of 
force-free magnetic field vectors is generally not force-free, but we will prove in 
the following (Equations (46) and (47)) that the sum of NLFFF solutions of the 
form of Equations (34)- (37)), which are force-free to second-order accuracy in a 
(or, more strictly, in [br sind]), have the property that their sum is also force- free 
to second-order in a. 

Let us first consider the condition of divergence-freeness. Since the divergence 
operator is linear, the superposition of a number of divergence-free fields is 
divergence-free also, 

V-B = V.(^B,) = ^(V.B,)==0. (45) 
i i 

While the divergence-free condition is exactly fulfilled for a potential field so- 
lution (Equation (4)), our quasi- forcefree approximation (Equations (34)- (37)) 
matches this requirement to second order in a, as the insertion of the solutions 
(Equations (34)-(37)) into the divergence expression (Equation (24)) shows. For 
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a quantitative measure of this level of accuracy we can also check numerical tests 
of the figure of merit (Section 3.3). 

Now, let us consider the condition of force-freeness. A force-free field has to 
satisfy Maxwell's equation (Equation (15)). Since we parameterized both the 
potential field and the non-potential field with a linear sum of iVm magnetic 
charges, the requirement would be, 

V X B = V X g B, = ^(V, X B,) = ^ a,{r)B, = a(r)B . (46) 

Generally, these three equations of the vector V x B cannot be fulfilled with a 
scalar function a{r) for a sum of force-free field components, unless the magnetic 
field volume consists of spatially separated force-free subvolumes. However, we 
can show the validity of the force-freeness equation (Equation (46)) to second- 
order accuracy in a. Note, that the nonlinear force-free parameter a is propor- 
tional to b (Equations (37) and (43)), which is defined in Equation (17), and 
thus we set second-order accuracy in b equal to second-order accuracy in a. The 
argument goes as follows. If we use spherical coordinates, the NLFFF solution 
of the radial component is of zeroth order, Br{r,6) oc 0(a°) (Equation (34)), 
the azimuthal component is of first order, B^{r,9) oc O(a^) (Equation (35)), 
and the neglected third component magnetic field component is of second-order, 
Be{r,e) (X O(a^) (as it can be shown by inserting Br and B^p into Equation (26). 
The curl of the magnetic field (Equations (25)-(27)) is then of first order for the 
radial component, [V x B],. oc aB,- oc 0{a^) (Equation (25)), to second order 
for the azimuthal component, [V x B]^ oc (aB^) oc 0{a^) (Equation (27)), and 
the remaining third component is of third-order, [V x 'B]g oc {aBg) oc 0{a^) 
(Equation (26)). 

Therefore, if we neglect second-order and higher-order terms, the divergence- 
free condition (Equation (46)), which generally has three equations for the three 
curl components, e.g., [V x Br]r, [V x B^]^, [V x Bg]g, reduces to one single 
equation for the radial component, [V x Br]r, which can be fulfilled with a scalar 
function a{r), 

^ [V X B],. _ [VxEfriB,]. _ E.t'iVxB,], _ E^=i»,B, 

Thus, we expect that the force-freeness is fulfilled to second-order accuracy 
O(a^) (or strictly speaking 0(&^)). We will demonstrate the near force-freeness 
of simulated examples in the next section. 

3. Simulations and Tests 

We are now going to simulate examples of the analytical nonlinear force-free 
solutions in order to visualize the magnetic topology and to quantify the accuracy 
of the divergence-free and force-free conditions. 
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Figure 6. Simulations of three line-of-sight magnetograms (left) and magnetic field lines 
projected into the x — y plane (left) and into the vertical x — z plane (right). Thre three 
cases include; (A) a single positive magnetic charge (first row), (B) a dipole produced by two 
magnetic charges with opposite polarity (second row), and (C) a quadrupole configuration 
(third row). See parameters in Table 1. Only field lines with magnetic fields above a 50% 
threshold of the maximum field strength are shown. 
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Figure 7. Simulations of three line-of-sight magnotograms (loft) and magnetic field lines of a 
non-potential model with currents are shown, projected into the x — y plane (left) and into the 
vertical x — z plane (right). The parameters of the three cases (D), (E), and (F) are identical 
to thocs of (A), (B), and (C), except for the addition of electric currents. 
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Figure 8. Simulations of tliree line-of-sight magnctograms (left) and magnetic field lines of 
a non-potential model with currents are shown, projected into the x — y plane (left) and into 
the vertical x — z plane (right). The three cases (G), (H), and (I) have each A^ni = 10 magnetic 
charges, with randomly chosen field strengths, locations, and electric currents. 
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3.1. Numerical Examples 

The simplest case is a single magnetic charge j ~ 1, which we illustrate as 
case A in Figure 6 (top row). We choose the following parameters: a magnetic 
field strength of Bi = 1000 G (gauss) at the solar surface directly above the 
buried charge, the location {xi,yi,Zi) = (0.1,0.0,0.95) for the buried charge, 
and a number of zero twist 6i = for the potential field case. We show the 
simulated line-of-sight magnetogram Bz{x, y) in Figure 1 (top left), which mimics 
an isolated sunspot. The pixel size of the magnetogram and the stepping size 
in the extrapolation along a field line is As = 0.004 solar radii (2800 km « 4", 
corresponding to the pixel size of SoHO/MDI magnetograms) . We extrapolate 
the field lines for every pixel that has a footpoint magnetic field strength above 
a threshold of 50% (> 500 G). The field lines point in radial direction away from 
the center of the buried magnetic charge, as it is expected for the potential field 
of an isolated sunspot (and defined in Equation (1)). 

The next basic example is a magnetic dipole, which can be represented in our 
model by a superposition of a pair of two magnetic charges with opposite polarity, 
as sketched in Figure 1. The case B shown in Figure 6 is simulated with equal, 
but oppositely signed magnetic field strengths [Bi = 1000 G, B2 — —1000 G) at 
mirrored positions (xi = 0.1, X2 = —0.1), otherwise we used the same parameters 
as in case A (yi = y2 = 0.0, ri = r2 = 0.95,^2 = 61 = 0.0). The magnetic field 
lines mimic the familiar structure of a dipole, which is parameterized here with 
8 free parameters (in the potential case). 

A quadrupolar configuration is simulated in case C (Figure 6, bottom), with 
translational symmetry {xi — 0.1, 2:2 = 0.05, X3 = —0.05, 2:4 = — 0.1;yi = 
0.1,2/2 = 0.05,2/3 ~ 0.1,1/4 ~ 0.05), equal depths (ri = 7-2 = ''3 = ''4 = 0.95), 
and ahernating field strengths (Bi = B3 = 1000,^2 = B4 = -1000 G). The 
quadrupolar configuration shows essentially two bipoles, each one with field lines 
that mostly connect within the same dipole domain, but a few intermediate field 
lines actually connect from one to the other domain. 

In Figure 7 we show the same three configurations as for the potential field 
model (A, B, C of Figure 6), but add electric currents caused by twisting, cor- 
responding to A'^twist = —0.5 turns for the single charge (case D) or first dipole 
(case E), and A^twist = 1.0 for the second dipole (case F), defined for a loop 
length of L = O.Itt solar radii. These amounts of twist correspond to force- free 
a-parameters of a = 2nNtv,ist/ L = —10 and —20 solar radius"^ {i.e., a = —1.43 
and —2.86 x 10~^° cm^^). Comparing the potential (Figure 6) and non-potential 
cases (Figure 7) shows clearly the differences that result from the presence of 
electric currents. The force-free field lines of a sunspot become distorted into 
spiral shapes (case D), the straight dipole becomes distorted into a sigmoid 
shape (case E) , and the quadrupolar configuration becomes also more distorted 
with sigmoid-like structures (case F). 

In Figure 8 we show a few more complicated cases (G, H, and I), consisting 
of A'ni = 10 magnetic charges, with random values chosen in the magnetic field 
range -1000 G < Bj < -1-1000 G, in positions -0.15 < Xj < 0.15 solar radii, 
—0.15 < yj < 0.15 solar radii, 0.95 < rj < 0.97 solar radii, and random twist in 
the range —3 < A'twist < +3 per L = O.Itt solar radii. The field lines displayed 
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Figure 9. Maps of the magnetic field components Bx{x,y), By(x,y), Bz{x,y) (left panels), 
the electric current density jz(x,y), and the force-free o-parameter (right panels). 



in Figure 8 demonstrate that a rich variety of signioid-shaped dipoles and inter- 
connecting niuhi-pole configurations can be generated with our quasi-force-free 
solutions, which mimic reahstic active regions observed in the solar corona. 
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Table 1. Figures of merit for nine simulations of nonlinear force-free field solutions, 
detailing the size of the 3D datacube, the number of magnetic charges (Afm), poten- 
tial or non-potential model (P and NP), the number of computed field lines TVf, the 
divergence- freeness Lj, the force- freeness Lf , and the computation times tcPU- 

Case Data cube Magnetic Field Divergence- Force- Computation 
charges lines freeness freeness time 

A^m Ni Ld Li tcpu (s) 



A 


51 


X 


51 


X 37 


1 (P) 


87 


0.0004 


0.0007 


0.078 


B 


51 


X 


51 


X 37 


2 (P) 


160 


0.0009 


0.0014 


0.309 


C 


51 


X 


51 


X 37 


4(P) 


159 


0.0015 


0.0019 


0.351 


D 


51 


X 


51 


X 37 


1 (NP) 


87 


0.0006 


0.0009 


0.083 


E 


51 


X 


51 


X 37 


2 (NP) 


160 


0.0007 


0.0010 


0.314 


F 


51 


X 


51 


X 37 


4 (NP) 


159 


0.0015 


0.0024 


0.414 


G 


51 


X 


51 


X 37 


10 (NP) 


336 


0.0012 


0.0058 


2.462 


H 


51 


X 


51 


X 37 


10 (NP) 


302 


0.0010 


0.0099 


1.764 


I 


51 


X 


51 


X 37 


10 (NP) 


217 


0.0018 


0.0133 


1.370 



3.2. Force-Free a-Parameter and Electric Current Maps 



In Figure 9 we show examples of various maps that can be generated to visualize 
a 3D vector field solution, for the case F of a quadrupolar configuration with 
currents. We show the following quantities in the image plane {x,y, z = 1 + As), 
which corresponds to an image plane near the solar surface: The three magnetic 
field vector component maps Bx{x,y), By{x,y), Bz{x,y) (Figure 9, left panels, 
the vertical electric current map jz{x, y) (Figure 9, top right panel), the nonlinear 
a-parameter a{x,y) (Figure 9, middle right panel), and the LOS magnetogram 
Bz{x,y) together with extrapolated field lines (Figure 9, bottom right panel). 
The Bz map shows most clearly the locations of the four buried magnetic charges 
that form two dipolar or a quadrupolar configuration. The magnetic polarization 
is also reflected in the and a-map. The B^ and the a-map show also the 
location of the neutral line, where numerical effects due to the limited spatial 
resolution become visible. 



3.3. Figures of Merit 



The degree of convergence towards a divergence-free magnetic field model solu- 
tion can be quantified by a measure that compares the average divergence V • B, 
which should be close to zero, with the gradient B/Ax of the magnetic field over 
a reference length scale Aa;, for instance a pixel of the computational grid. The 
average deviation can then be defined by (see also Wheatland et ai, (2000) or 
Schrijver et ai, (2006)), 
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Similarly, the force-freeness can be quantified by the ratio of the Lorentz force, 
(j X B) = (V X B) X B to the normalization constant B^/Ax, 



where B ^ \B\. 

We calculated these figure of merit quantities for the nine cases simulated in 
Figures 6-9. The values are listed in each of the panels in Figures 6-8 and listed 
in Table I. The potential-field cases (A, B, and C) are found to have a figure 
of merit in the range of Ld = 0.0009 ± 0.0006 for the divergence-freeness, and 
L{ = 0.0014 ± 0.0006 for the forcc-frccness. The non-potential field cases (D, E, 
F, G, H, and I) have values in similar ranges of = 0.0009 ± 0.0005 for the 
divergence-freeness, and L{ = 0.0100 ± 0.0080 for the force-freeness. We find no 
tendency that this figure of merit depends on the number of magnetic charges 
or some other model parameters. The fact that our quasi-force free analytical 
solutions perform equally well as standard NLFFF codes described in Schrijver 
eq al. (2006) tells us that the inaccuracy of the analytical approximation is 
commensurable or even smaller than the numerical uncertainty of other NLFFF 
codes. However, since our analytical solution provides an explicit formulation 
of nonlinear force-free fields, it can be computed much faster than the standard 
NLFFF codes, and still provides approximate solutions with acceptable accuracy 
(to second order). The computation time of the analytical solutions for the cases 
shown in Figures 6-8 amounts to about 1 s (on a recent Mac computer: Mac 
OS X, 2 X 3.2 GHz Quad-Core Intel Xeon, Memory 32 GB 800 MHz DDR2 FB- 
DIMM), while standard iterative NLFFF codes need several hours to converge 
to a single NLFFF solution. 

4. Discussion and Conclusions 

The coronal magnetic field has generally been computed by extrapolation from 
lower boundary data in form of photospheric magnetograms Bz{x,y,z = Zph) 
or vector-magnetograph data B(a;, y), using a numerical extrapolation algo- 
rithm that fulfills the conditions of force-freeness (V-B) and divergence-freeness 
V X B = a(r)B, where a(r) is a scalar function in space r. These extrapolation 
algorithms are very computing-intensive, because a good solution requires many 
iterations on a large computational 3D-grid that has sufficient spatial resolution 
to resolve the relevant magnetic field gradients. The accuracy of these numerical 
solutions depends very much on the noise in boundary vector magnetic field 
data as well as on deviations of photospheric fields from a force-free state. Re- 
cent stereoscopic triangulation of coronal loops has demonstrated a considerable 
mismatch between the extrapolated fields and the actual coronal loops, which 
cannot easily be reconciled with extrapolation algorithms, since they have only a 
very limited degree of freedom within the noise of the boundary data. Moreover, 
since each NLFFF solution is very time-consuming to compute, these algorithms 
arc not suitable for forward-fitting. 




(49) 



SOLA: ms.tex; 13 July 2012; 0:13; p. 20 



Nonlinear Force-Free Magnetic Field 



The forward-fitting of magnetic field solutions to observed data requires a 
faster algorithm to compute many NLFFF solutions for variable boundary data 
or for coronal constraints as given by stereoscopic 3D reconstructions. The fastest 
computational way would be an explicit analytical solution for the coronal field 
vectors B(r) as a function of some suitable parameterization of the boundary 
data or coronal constraints. There exist some analytical solutions of nonlinear 
forcc-frcc fields, such as a class of solutions in terms of Lcgendre polynomi- 
als (Low and Lou, 1990), which is characterized by some spatial symmetry 
and has been used to test numerical extrapolation algorithms {e.g., DeRosa 
et al, 2009; Malanushenko et ai, 2009). However, to our knowledge, the class 
of analytical NLFFF solutions of Low and Lou (1990) has never been applied 
to forward-fitting of observed data, such as line-of-sight magnetograms, vector 
magnetograph 3D data, or to stcrcoscopically triangulated loops. Moreover, the 
special class of NLFFF solutions derived in Low and Lou (1990) correspond to 
harmonics of Lcgendre polynomials, which have a high degree of symmetry that 
does not match realistic observations of active regions, and thus is not suitable 
for forward-fitting to real data. 

What we need to model observed solar magnetic data with high accuracy is: 
(1) an explicit formulation of an analytical NLFFF solution; (2) a parameteriza- 
tion of the NLFFF solution with a sufficient large number of free parameters that 
can be forward-fitted to data and converges close to observations; and (3) a fast 
computation algorithm that can perform many interations without computing- 
intensive techniques. Hence, such a project consists of developing a suitable 
analytical formulation first, and then to implement the analytical solutions into 
a forward-fitting code. In this paper we have undertaken the first step. We started 
with a potential-field parameterization in terms of 7V,„ buried magnetic charges, 
which is defined by 4iVm free parameters that can easily be extracted from an 
observed line-of-sight magnetogram Bz{x, y) with arbitrary accuracy, as demon- 
strated in two recent studies (Aschwanden and Sandman, 2010; Aschwanden 
et ai, 2012a). The key concept of this potential-field representation is that an 
arbitrary complex 3D magnetic field can be decomposed into a finite number 
of elementary magnetic field components, where each one simply consists of a 
quadratically decreasing radial field of a buried magnetic charge. Divergence- 
freeness is conserved due to the linearity in the superposition of elementary 
field components. In a next step we extended the elementary potential-field 
component to a nonpotential-field component by adding a uniform twist that 
can be parameterized by the force-free a-parameter. Such an elementary non- 
potential field component requires five free parameters, consisting of the four 
potential-field parameters plus the force-free a-parameter. We derived an explicit 
analytical formulation of the radial Br(r,9) and azimuthal field vector B^{r,9) 
that represents an approximative solution of the divergence-free and force-free 
condition to second order (oc a^). This solution is very accurate for weakly non- 
potential fields and converges to the potential field solution for a = 0. In analogy 
to the potential-field representation, we represent a general non-potential field 
solution with a superposition of elementary non-potential field components and 
prove that the divergence- frecness and force- freeness is conserved to second-order 
accuracy in our NLFFF approximation. 
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We calculated some examples of potential and non-potential fields that mimic 
an isolated sunspot, a dipolar and a quadrupolar configuration, as well as more 
complex multi-polar configurations. The examples show that the magnetic field 
of arbitrary complex active regions can be represented with our parameteriza- 
tion. Increasing the force- free a-parameter distorts circular field lines into helical 
and sigmoid-shapcd geometries. Our parameterization allows one to compute 
either field lines (starting from arbitrary locations), 3D datacubcs of magnetic 
field vectors, of maps of the force-free a-parameter and electric current jz (Figure 
9). We tested the figures of merit for divergence-freeness and force- freeness, 
which amount to Ld $ 10~^ and Lf < 10~^. The examples demonstrate also the 
computing speed of this algorithm, which amounts to the order of w 1 s for 
a computation grid that encompasses a typical active region with the spatial 
resolution of MDI. Thus, we envision that a full-fietched forward-fitting code 
can converge within a few seconds to a few minutes, depending on the number 
of iterations and number of magnetic field components. 

Where do we go from here? The next step is the development of a forward- 
fitting code that uses the magnetic field parameterization described here (see 
Paper II). We envision the applications to at least three different sets of con- 
straints, requiring three different versions of forward-fitting codes: (i) line-of- 
sight magnetograms Bz{x, y) and 3D coordinates [x{s), y{s), z{s)] of stereoscop- 
ically triangulated loops; (ii) line-of-sight magnetograms Bz{x,y) and 2D coor- 
dinates [(a:(s), y{s)] of traced loops; and (ii) vcctor-magnetograph data [Bx{x, y), 
By{x,y), Bz{x,y)]. The first application requires STEREO data, while the sec- 
ond one can be obtained from any EUV imager {e.g., AIA/SDO, TRACE, 
EIT/ SOHO). The third application can be conducted with the new HMI/SDO 
data and is equivalent to other NLFFF extrapolation codes without coronal 
constraints, while the first two use coronal tracers and alleviate the force-free 
assumption of photospheric data. We envision that these three applications will 
reveal insights into a number of crucial questions in a novel way. 

There is a large number of physical problems and issues that can be addressed 
with the anticipated forward-fitting code, such as: (i) The force-frecness of the 
photosphere; (ii) the accuracy of NLFFF solutions; (iii) the spatial distribution 
of electric currents in active regions; (iv) the temporal evolution of currents 
before and during fiares; (v) the spatial distribution of current dissipation and 
coronal heating; (vi) helicity injection; (vii) the 3D geometry of coronal loops 
which is needed for hydrodynamic modeling; (viii) scaling laws of the volumetric 
heating function with other physical parameters; (ix) tests of the magnetic field 
strength inferred from coronal seismology, eic. There exists hardly a phenomenon 
in the solar corona that can be modeled without the knowledge of the coronal 
magnetic field. 

Appendix A: The Gold-Hoyle Flux Rope 

A simple geometry of a force-free field structure is the Gold-Hoyle flux rope 
(Gold and Hoyle, 1960), which consists of a curved axis with helical field lines 
curved around the axis (Figure 10). While the stretched version of a flux rope 
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Figure 10. Cartoon of Gold-Hoylc flux rope. 

with a straight twist axis has the exact force-free solution of a uniformly twisted 
flux tube (Section 2.2), the curved version of the Gold-Hoyle flux rope is subject 
to curvature forces due to the gradient of the magnetic field across the flux rope 
diameter and has a modified force-free solution. 

In order to explore the limitations of our force-free field parameterization we 
attempt here to model such a Gold-Hoyle flux rope. We use the coordinates 
(xqjO, zo) and (— xo,0, zq) with xq = 0.1 and zq = 0.985 solar radii (marked 
with diamonds in Figure 11) and extrapolate field lines B{s) with our method, 
starting from the apex position (0, 0, Za) with Za = 1.1, for a set of six cases with 
various force-free parameters ai — q;2, where the a's associated with the twist 
axis of each buried charge are defined by a = 27r7Vtwist/-^, with the loop length 
L = 27rxo = 0.314 and the number of twist turns iVtwist = 0, 1, 5 (indicated 
with N = 0, 5 in Figure 11). The case ~ corresponds to the potential 
field case, yielding a coplanar elliptical loop shape. The case iV = 1 represents 
a slightly twisted field line that has a sigmoid shape and is a quasi-force-free 
solution. The cases with = 2, 5 are strongly twisted field lines and may be 
less force-free, since the neglected terms could be significant. 

Obviously we cannot reproduce the exact shape of the Gold-Hoylc flux rope 
as shown in Figure 10 (with about seven twist turns) with our choice of parame- 
terization. The reason lies in the geometric constraints of the twist axis, which is 
semi-circular in the case of the Gold-Hoyle model, but consists of vertical twist 
axes in our parameterization. So, this counter-example clearly demonstrates the 
limitations of our parameterization. Nevertheless, although the cartoon with the 
Gold-Hoyle geometry is very popular, especially for interplanetary flux ropes 
and CMEs, it is not clear whether such Gold-Hoyle type geometries are found in 
loops in the lower corona, and whether the Gold-Hoyle geometry corresponds to 
an exact force- free solution. It is conceivable that the Sun exerts rotational stress 
mostly in the photosphere {i.e., rotating sunspots), which propagates in vertical 
direction along the field lines, but does not necessarily lead to a uniformly twisted 
circular flux tube as shown in Figure 10, because the magnetic field drops rapidly 
with with height (for magnetic charges with small sub-photosphcric depths), 
and thus the magnetic stress is not uniformly distributed along a semi-circular 
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Figure 11. Dipolar field lines with various numbers of twisting turns: N = (potential field 
line), stable sigmoid (A'' = 1; solid line), and unstable sigmoids [N = 2,..., 5; dashed lines), 
aecording to our parameterization of point charges with twisted vertical axes. Note that the 
limit of large twist numbers does not turn into a Gold-Hoyle flux rope (Figure 10) with our 
parameterization. 



potential field line as envisioned in the Gold-Hoyle scenario. However, for a case 
with a near-constant magnetic field strength B{s) along a potential field line, 
we would expect a uniform twist as outlined in the Gold-Hoyle case. 

On the other side, strongly twisted fiux tubes with a twist larger than about 
1.25 full turns are unstable due to the kink instability and may erupt, which 
is another reason why multiply twisted fiux tubes are unlikely to be found in 
active regions. Even Gold and Hoyle (1960) found a critical twist number of 
*&twist ~ 2.497r (A^twist = <i'twist/27r < 1.25) above which no equilibrium exists, 
which is also confirmed by recent MHD simulations {e.g., Torok and Kliem, 
2003). Thus, the Gold and Hoyle flux rope case may not be relevant for modeling 
magnetic fields in stable active regions. Nevertheless, more general parameteriza- 
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tions could be anticipated in future work, such as twist axes that foUow potential 
field lines, rather than vertical used in our parameterization to minimize 

the number of free parameters. 
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